

%% data

dChar = readtable('../../../raw_data/district_characteristics.csv');
e09 = readtable('../../../raw_data/election_results_2009.csv');
e12 = readtable('../../../raw_data/election_results_2012.csv');
cs12 = readtable('../../../raw_data/campaign_spending_2012.csv');

% order observations by candidate-coalition choice
d0 = find(e12.candCM==0); N0 = size(d0,1);
d1 = find(e12.candCM==1); N1 = size(d1,1);
d2 = find(e12.candCM==2); N2 = size(d2,1); N = N0+N1+N2;
dChar = dChar([d0;d1;d2],:);
e09 = e09([d0;d1;d2],:); e12 = e12([d0;d1;d2],:);
cs12 = cs12([d0;d1;d2],:);
clear d0 d1 d2

% district indicator/aggregation matrix
A = [1:N,1:N,1:N0,N0+N1+(1:N2),1:N0+N1,1:N]';
A = 1 * (A==A');

% lagged vote shares
P = [e09.MP, ...
    e09.NA, ...
    e09.PVEM + 0.5 * e09.PM, ...
    e09.PRI + 0.5 * e09.PM, ...
    e09.PAN] ./ e09.registered;

% 2012 vote shares
S = [e12.MP ./ e12.registered; ...
    e12.NA ./ e12.registered; ...
    e12.PVEM(1:N0,1) ./ e12.registered(1:N0,1); ...
    (e12.PRI(N0+N1+(1:N2),1) + e12.PVEM(N0+N1+(1:N2),1) + e12.CM(N0+N1+(1:N2),1)) ./ e12.registered(N0+N1+(1:N2),1); ...
    e12.PRI(1:N0,1) ./ e12.registered(1:N0,1); ...
    (e12.PRI(N0+(1:N1),1) + e12.PVEM(N0+(1:N1),1) + e12.CM(N0+(1:N1),1)) ./ e12.registered(N0+(1:N1),1); ...
    e12.PAN ./ e12.registered];


%% baseline specification

% demographics
D = [dChar.femaleHH,dChar.over60,dChar.rural];
dnames = {'female','over60','rural'};
KD = size(D,2);

% menu fixed effects and demographics
X = [blkdiag(ones(N0,1),ones(N1,1),ones(N2,1)),D];
X = blkdiag(X,X,[blkdiag(ones(N0,1),ones(N2,1)),[D(1:N0,:);D(N0+N1+(1:N2),:)]],...
    [blkdiag(ones(N0,1),ones(N1,1)),D(1:N0+N1,:)],X);
NN = size(X,1);

% add log-lagged vote shares
PP = zeros(NN,3);
PP([1:N0,N+(1:N0),2*N+(1:N0),2*N+N0+N2+(1:N0),3*N+N0+(1:N0)],1) = reshape(log(P(1:N0,:)),5*N0,1);
PP([N0+(1:N1),N+N0+(1:N1),2*N+N0+N2+N0+(1:N1),3*N+N0+N0+(1:N1)],2) = reshape(log([P(N0+(1:N1),1:2),sum(P(N0+(1:N1),3:4),2),P(N0+(1:N1),5)]),4*N1,1);
PP([N0+N1+(1:N2),N+N0+N1+(1:N2),2*N+N0+(1:N2),3*N+N0+N0+N1+(1:N2)],3) = reshape(log([P(N0+N1+(1:N2),1:2),sum(P(N0+N1+(1:N2),3:4),2),P(N0+N1+(1:N2),5)]),4*N2,1);
X = [X,sum(PP,2)];

% add governors
G = zeros(NN,1);
G(1:N,:) = strcmp(dChar.gov12,'MC') + strcmp(dChar.gov12,'PRD');
G(N+(1:N),:) = strcmp(dChar.gov12,'NA');
G(2*N+(1:N0),:) = strcmp(dChar.gov12(1:N0),'PVEM');
G(2*N+N0+(1:N2),:) = strcmp(dChar.gov12(N0+N1+(1:N2)),'PVEM') + strcmp(dChar.gov12(N0+N1+(1:N2)),'PRI');
G(2*N+N0+N2+(1:N0),:) = strcmp(dChar.gov12(1:N0),'PRI');
G(2*N+2*N0+N2+(1:N1),:) = strcmp(dChar.gov12(N0+(1:N1)),'PVEM') + strcmp(dChar.gov12(N0+(1:N1)),'PRI');
G(3*N+N0+(1:N),:) = strcmp(dChar.gov12,'PAN');
X = [X,G];
K0 = size(X,2);

% (endogenous) campaign spending
C = [cs12.MP;...
    cs12.NA;...
    cs12.PVEM(1:N0,1);cs12.CM(N0+N1+(1:N2),1);...
    cs12.PRI(1:N0,1);cs12.CM(N0+(1:N1),1);...
    cs12.PAN];
XC = [X,C,C.^2];
K = size(XC,2);


%% baseline specification: voting second stage

X2S = [blkdiag(ones(N1,1),ones(N2,1)),D(N0+1:end,:)];
X2S = blkdiag(X2S,X2S);
X2S=[X2S,[log(P(N0+1:end,3));log(P(N0+1:end,4))],...
    [strcmp(dChar.gov12(N0+1:end),'PVEM');strcmp(dChar.gov12(N0+1:end),'PRI')]];

% OLS
Y2S = [log(e12.PVEM(N0+1:end)) - log(e12.CM(N0+1:end));...
    log(e12.PRI(N0+1:end)) - log(e12.CM(N0+1:end))];
B2S = X2S \ Y2S;
xi2S = Y2S - X2S * B2S;
VB2S = (X2S' * X2S) \ ...
    (X2S' * diag(xi2S.*xi2S) * X2S) / (X2S' * X2S); % heteroskedasticity-robust variance
seB2S = sqrt(diag(VB2S)); % standard errors
rB2S = [B2S,seB2S,2*(1-normcdf(abs(B2S./seB2S)))];

disp(' ')
disp('for Table D2 (column (II)):')
disp(print_results_govs(round(rB2S,3),dnames,[],[],1))


%% baseline specification: BLP

load('../../voting_stage.mat','xi0','rBblp','Z0')
KZ = size(Z0,2);

% SGI nodes and weights
[nu,wnu] = nwspgr('KPN',2,6);

% BLP
W = Z0' * diag(xi0.*xi0) * Z0; % GMM weighting matrix
Bblp = [rBblp(1:K0-1,1);0;rBblp(K0:end,1)];
JPattern = [ones(length(S),length(Bblp)),A,sparse(length(S),KZ)];
HPattern = sparse(blkdiag([ones(length(Bblp),length(Bblp)+length(S));ones(length(S),length(Bblp)),A],ones(KZ)));
optVS = optimset('Display','off','JacobPattern',JPattern,'HessPattern',HPattern,...
    'HessFcn',@(BB,lambda)HessL(BB,lambda,XC,KD,A,[N0,N1,N2],W,nu,wnu));
[B,~,exit] = knitro_nlp(@(BB)Q_BLP(BB,length(Bblp),length(S),W),[Bblp;xi0;Z0'*xi0],[],[],...
    [sparse(KZ,length(Bblp)),Z0',-speye(KZ)],sparse(KZ,1),...
    [-Inf(length(Bblp)-2,1);zeros(2,1);-Inf(length(S)+KZ,1)],Inf(length(Bblp)+length(S)+KZ,1),...
    @(BB)Shares(BB,S,XC,KD,A,[N0,N1,N2],nu,wnu),[],optVS,'knitro.opt');
xi = B(length(Bblp)+(1:length(S)),1);

% covariance matrix
[~,~,~,Ghat] = Shares(B,S,XC,KD,A,[N0,N1,N2],nu,wnu);
Ghat = Ghat';
Ghat = - [XC,Ghat(:,K+2+(1:length(S)))\Ghat(:,K+(1:2))];
Ghat = Z0' * Ghat;
VBblp = inv(Ghat'*((Z0'*diag(xi.*xi)*Z0)\Ghat)); % robust variance
seB1 = sqrt(diag(VBblp)); % standard errors
Bblp = B(1:length(Bblp),1);
rBblp = [Bblp,seB1,2*(1-normcdf(abs(Bblp./seB1)))];

disp(' ')
disp('for Table D1 (column (II)):')
disp(print_results_govs(round(rBblp,3),dnames,[],1,0))


%% region FE

% demographics
D = [dChar.region<=2,dChar.region==3,dChar.femaleHH,dChar.over60,dChar.rural];
rfe = {'north','southeast'};
KD = size(D,2);

% menu fixed effects and demographics
X = [blkdiag(ones(N0,1),ones(N1,1),ones(N2,1)),D];
X = blkdiag(X,X,[blkdiag(ones(N0,1),ones(N2,1)),[D(1:N0,:);D(N0+N1+(1:N2),:)]],...
    [blkdiag(ones(N0,1),ones(N1,1)),D(1:N0+N1,:)],X);

% add log-lagged vote shares
X = [X,sum(PP,2)];

% add governors
X = [X,G];
K0 = size(X,2);

% (endogenous) campaign spending
XC = [X,C,C.^2];
K = size(XC,2);


%% region FE: voting second stage

X2S = [blkdiag(ones(N1,1),ones(N2,1)),D(N0+1:end,:)];
X2S = blkdiag(X2S,X2S);
X2S=[X2S,[log(P(N0+1:end,3));log(P(N0+1:end,4))],...
    [strcmp(dChar.gov12(N0+1:end),'PVEM');strcmp(dChar.gov12(N0+1:end),'PRI')]];

% OLS
B2S = X2S \ Y2S;
xi2S = Y2S - X2S * B2S;
VB2S = (X2S' * X2S) \ (X2S' * diag(xi2S.*xi2S) * X2S) / (X2S' * X2S); % heteroskedasticity-robust variance
seB2S = sqrt(diag(VB2S)); % standard errors
rB2S_RFE = [B2S,seB2S,2*(1-normcdf(abs(B2S./seB2S)))];

disp(' ')
disp('for Table D2 (column (VI)):')
disp(print_results_govs(round(rB2S_RFE,3),dnames,rfe,[],1))


%% region FE: BLP

load('../../voting_stage.mat','xi0rfe','rBblp_RFE','Z0rfe')
KZ = size(Z0rfe,2);

% BLP
W = Z0rfe' * diag(xi0rfe.*xi0rfe) * Z0rfe; % GMM weighting matrix
Bblp = [rBblp_RFE(1:K0-1,1);0;rBblp_RFE(K0:end,1)];
JPattern = [ones(length(S),length(Bblp)),A,sparse(length(S),KZ)];
HPattern = sparse(blkdiag([ones(length(Bblp),length(Bblp)+length(S));ones(length(S),length(Bblp)),A],ones(KZ)));
optVS = optimset('Display','off','JacobPattern',JPattern,'HessPattern',HPattern,...
    'HessFcn',@(BB,lambda)HessL(BB,lambda,XC,KD,A,[N0,N1,N2],W,nu,wnu));
[B,~,exit_RFE] = knitro_nlp(@(BB)Q_BLP(BB,length(Bblp),length(S),W),[Bblp;xi0rfe;Z0rfe'*xi0rfe],[],[],...
    [sparse(KZ,length(Bblp)),Z0rfe',-speye(KZ)],sparse(KZ,1),...
    [-Inf(length(Bblp)-2,1);zeros(2,1);-Inf(length(S)+KZ,1)],Inf(length(Bblp)+length(S)+KZ,1),...
    @(BB)Shares(BB,S,XC,KD,A,[N0,N1,N2],nu,wnu),[],optVS,'knitro.opt');
xi = B(length(Bblp)+(1:length(S)),1);

% covariance matrix
[~,~,~,Ghat] = Shares(B,S,XC,KD,A,[N0,N1,N2],nu,wnu);
Ghat = Ghat';
Ghat = - [XC,Ghat(:,K+2+(1:length(S)))\Ghat(:,K+(1:2))];
Ghat = Z0rfe' * Ghat;
VBblp = inv(Ghat'*((Z0rfe'*diag(xi.*xi)*Z0rfe)\Ghat)); % robust variance
seB1 = sqrt(diag(VBblp)); % standard errors
Bblp = B(1:length(Bblp),1);
rBblp_RFE = [Bblp,seB1,2*(1-normcdf(abs(Bblp./seB1)))];

disp(' ')
disp('for Table D1 (column (VII)):')
disp(print_results_govs(round(rBblp_RFE,3),dnames,rfe,1,0))


%%

save('voting_stage_govs.mat')



